******************************************************************
** Zhirnov, Moral, Sedashov - Taking Distributions Seriously *****
** Replication Do File: Figures (January 20, 2022) **********
******************************************************************
clear all
set more off
set matsize 800
set seed 123456789

set scheme s1mono
*** set scheme mmoral5 /* The scheme file is available from the authors */

timer clear
timer on 1
 
ssc install moremata, replace
ssc install estout, replace
net install grc1leg2, from("http://digital.cgdev.org/doc/stata/MO/Misc") replace force

timer on 2

******************************************************************
** Functions *****************************************************
******************************************************************

** Marginal Effects
mata
real matrix me_byrow(coef, X, Z) {
	dydx=me(coef, X, Z)
	means=mean(dydx)'
	ra=mm_quantile(dydx, 1, (0.025 \ 0.975))'						/* Confidence level can be changed here */
	return((means,ra))
}

real matrix me_wt(coef, X, Z, group_id, weight) {
	dydx=me(coef, X, Z)
	groups=uniqrows(group_id)
	wtm=J(cols(dydx), rows(groups), .)
	obs=J(rows(groups),1,.)
	for (i=1; i<=rows(groups); i++) {
		wtmc=(group_id:==groups[i]):*weight
		obs[i]=sum(wtmc)
		wtm[.,i]=wtmc/sum(wtmc)
		}
	dydxw=dydx*wtm
	means=mean(dydxw)'
	ra=mm_quantile(dydxw, 1, (0.025 \ 0.975))'	
	return((groups,obs,means,ra))
}
end

** Heatmap
capture program drop heatm
program heatm
args filename name1 name2 varname1 varname2 gform mform graph_name
use `filename', clear
gen significant=(lb>0 & ub>0)|(lb<0 & ub<0)
qui sum `varname1' [fw=obs]
loc x1mean: disp %9.`gform'f r(mean)
loc x1max: disp %9.`gform'f r(max)
loc x1min: disp %9.`gform'f r(min)
loc s1=round((`x1max'-`x1min')/4, `=1/10^`gform'')				/* Number of ticks on the y axis can be changed here or changing the `ylab' option below */		
qui sum `varname2' [fw=obs]
loc x2mean: disp %9.`gform'f r(mean)
loc x2max: disp %9.`gform'f r(max)
loc x2min: disp %9.`gform'f r(min)
loc s2=round((`x2max'-`x2min')/4, `=1/10^`gform'')				/* Number of ticks on the y axis can be changed here or changing the `xlab' option below */

qui sum dydx_dx1
loc dvmax: disp %9.`mform'f r(max)
loc dvmin: disp %9.`mform'f r(min)
loc step=round((`dvmax'-`dvmin')/3, `=1/10^`mform'')			/* Number of cutpoints can be changed here */

gen counter=_n
qui sum counter
loc extra1=`=r(max)'+1
loc extra2=`=r(max)'+2
loc extra3=`=r(max)'+3
loc extra4=`=r(max)'+4
set obs `extra4'

qui sum obs
replace significant=1 in `extra1'/`extra2'
replace significant=0 in `extra3'/`extra4'
replace obs=`=r(min)' in `extra1'/`extra4'
replace obs=`=r(max)' in `extra2'/`extra3'

if abs(`dvmax') > abs(`dvmin') {								/* The starting and end colors and intensity can be changed here */
	loc scolr="yellow*.25"
	loc ecolr="red*.95"
	}
else {
	loc scolr="red*.95"
	loc ecolr="yellow*.25"
}

twoway histogram `varname2' [fw=obs], frac ysca(alt reverse) ylab(#4, nogrid labsize(medsmall)) xlab(`x2min'(`s2')`x2max' `x2mean', grid gmax labsize(medsmall)) fysize(20) fcolor(black%95) lwidth(vthin) lcolor(white%25) xtitle(`name2', size(medsmall)) ytitle("") nodraw name(hy, replace)
twoway histogram `varname1' [fw=obs], frac xsca(alt reverse) horiz 	xlab(#4, nogrid labsize(medsmall)) ylab(`x1min'(`s1')`x1max' `x1mean', grid gmax labsize(medsmall)) fxsize(20) fcolor(black%95) lwidth(vthin) lcolor(white%25) ytitle(`name1', size(medsmall)) xtitle("") nodraw name(hx, replace)
twoway (contour dydx_dx1 `varname1' `varname2' if dydx_dx1!=., levels(100) crule(linear) scolor(`scolr') ecolor(`ecolr') zlab(`dvmin'(`step')`dvmax', labsize(medsmall))) /*
*/ (scatter `varname1' `varname2' [fw=obs] if significant==0, msymbol(oh) mlcolor(black%95) mlwidth(vthin) msize(*.25)) /*
*/ (scatter `varname1' `varname2' [fw=obs] if significant==1, msymbol(o) mfcolor(black%95) mlwidth(none) msize(*.25)), /*
*/ xsca(alt) ysca(alt) xtitle("") ytitle("") ztitle("") ylab(`x1min'(`s1')`x1max' `x1mean', grid gmax labsize(medsmall)) xlab(`x2min'(`s2')`x2max' `x2mean', labsize(medsmall) grid gmax) /* 
*/ legend(off) clegend(title("Effect Size", size(medsmall) pos(12) justification(right)) ring(0) width(5) height(25)) nodraw name(yx, replace)
gr combine hx yx hy, hole(3) imargin(zero) scale(1.1) xsize(5.5) ysize(5.5)
gr export `graph_name'.eps, replace
end

** Contour Plot
capture program drop contourp
program contourp
args filename name1 name2 varname1 varname2 gform mform graph_name
use `filename', clear
gen significant=(lb>0 & ub>0)|(lb<0 & ub<0)
qui sum `varname1' [fw=obs]
loc x1max: disp %9.`gform'f r(max)
loc x1min: disp %9.`gform'f r(min)
loc s1=round((`x1max'-`x1min')/4, `=1/10^`gform'')				/* Number of ticks on the y axis can be changed here or using the `ylab' option below */
qui sum `varname2' [fw=obs]
loc x2max: disp %9.`gform'f r(max)
loc x2min: disp %9.`gform'f r(min)
loc s2=round((`x2max'-`x2min')/4, `=1/10^`gform'')				/* Number of ticks on the y axis can be changed here or changing the `xlab' option below */

gen counter=_n
qui sum counter
loc extra1=`=r(max)'+1
loc extra2=`=r(max)'+2
loc extra3=`=r(max)'+3
loc extra4=`=r(max)'+4
set obs `extra4'

qui sum obs
replace significant=1 in `extra1'/`extra2'
replace significant=0 in `extra3'/`extra4'
replace obs=`=r(min)' in `extra1'/`extra4'
replace obs=`=r(max)' in `extra2'/`extra3'

qui sum dydx_dx1,detail
loc locut=`r(min)' +(`r(max)'-`r(min)')*1/4
loc medcut=`r(min)'+(`r(max)'-`r(min)')*2/4
loc hicut=`r(min)' +(`r(max)'-`r(min)')*3/4

if abs(`hicut') > abs(`locut') {
	loc colr="white*.5 yellow*.5 orange*.5 red*.5"	/* The starting and end colors, and intensity can be changed here */
	}
	else {
		loc colr="red*.5 orange*.5 yellow*.5 white*.5"
}
twoway (contour dydx_dx1 `varname1' `varname2' if dydx_dx1!=., ccuts(`locut' `medcut' `hicut') ccolors(`colr')) /*
*/ (scatter `varname1' `varname2' [fw=obs] if significant==0, msymbol(oh) mlcolor(black%95) mlwidth(vthin) msize(*.25)) /*
*/ (scatter `varname1' `varname2' [fw=obs] if significant==1, msymbol(o) mfcolor(black%95) mlwidth(none) msize(*.25)), /*
*/ xtitle(`name2', size(medsmall)) ytitle(`name1', size(medsmall)) ztitle("") ylab(`x1min'(`s1')`x1max', labsize(medsmall) grid gmax) xlab(`x2min'(`s2')`x2max', labsize(medsmall) grid gmax) ysize(5.5) xsize(5.5) /*
*/ zlabel(`r(min)' `locut' `medcut' `hicut' `r(max)', format(%9.`mform'f) labsize(medsmall)) legend(off) clegend(title("Effect Size", size(medsmall) pos(12) justification(right)) width(5) height(25)) name(yx, replace)
gr export `graph_name'.eps, replace
end

** DAME Plot
capture program drop dame_plot
program dame_plot
args filename name1 name2 varname2 v2alias mform graph_name
use `filename', clear
gen significant=(lb>0 & ub>0)|(lb<0 & ub<0)
gen counter=_n
qui sum counter
loc extra1=`=r(max)'+1
loc extra2=`=r(max)'+2
set obs `extra2'
qui sum obs
replace significant=0 in `extra1'/`extra2'
replace obs=`=r(min)' in `extra1'
replace obs=`=r(max)' in `extra2'
qui sum ub
loc dvmax: disp %9.`mform'f r(max)
qui sum lb
loc dvmin: disp %9.`mform'f r(min)
qui sum ubm
loc dvmaxm: disp %9.`mform'f r(max)
qui sum lbm
loc dvminm: disp %9.`mform'f r(min)
if `dvmax'>`dvmaxm' {
	loc total_max=`dvmax'
	}
else {
		loc total_max=`dvmaxm'
	}
if `dvmin'>`dvminm' {
		loc total_min=`dvminm'
	}
else {
		loc total_min=`dvmin'
}
loc step=round((`total_max'-`total_min')/4, `=1/10^`mform'')
qui sum significant if dydx_dx1!=.
if r(max)==0 {
	twoway (line dydx_mdx1 `varname2', lpattern(solid) lwidth(vthin) lcolor(black%50)) ///
	(rline lbm ubm `varname2', lpattern(dash) lwidth(vthin) lcolor(black%50)) ///
	(rspike lb ub `v2alias', lcolor(black%95) lwidth(thin)) ///
	(scatter dydx_dx1 `v2alias' if significant==0 [fw=obs], msymbol(oh) mlcolor(black%95) mlwidth(vthin) msize(*.25)), ///
	ytitle(`name1', height(4) size(medsmall)) xtitle(`name2', height(4) size(medsmall)) xlab(, labsize(medsmall)) ///
	ylab(`total_min'(`step')`total_max', labsize(medsmall) gmax) yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5) legend(off)
	gr export `graph_name'.eps, replace
}
if r(min)==1 {
	twoway (line dydx_mdx1 `varname2', lpattern(solid) lcolor(gs10%95)) ///
	(rline lbm ubm `varname2', lpattern(dash) lcolor(gs10%95)) ///
	(rspike lb ub `v2alias', lcolor(black%95) lwidth(thin)) ///
	(scatter dydx_dx1 `v2alias' if significant==1 [fw=obs], msymbol(o) mfcolor(black%95) mlwidth(none) msize(*.25)), ///
	ytitle(`name1', height(4) size(medsmall)) xtitle(`name2', height(4) size(medsmall)) xlab(, labsize(medsmall)) ///
	ylab(`total_min'(`step')`total_max', labsize(medsmall) gmax) yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5) legend(off)
	gr export `graph_name'.eps, replace
}
if r(min)==0 & r(max)==1 {
	twoway (line dydx_mdx1 `varname2', lpattern(solid) lwidth(vthin) lcolor(black%50)) ///
	(rline lbm ubm `varname2', lpattern(dash) lwidth(vthin) lcolor(black%50)) ///
	(rspike lb ub `v2alias', lcolor(black%95) lwidth(thin)) ///
	(scatter dydx_dx1 `v2alias' if significant==0 [fw=obs], msymbol(oh) mlcolor(black%95) mlwidth(vthin) msize(*.25)) ///
	(scatter dydx_dx1 `v2alias' if significant==1 [fw=obs], msymbol(o) mfcolor(black%95) mlwidth(none) msize(*.25)), /// 
	ytitle(`name1', height(4) size(medsmall)) xtitle(`name2', height(4) size(medsmall)) xlab(, labsize(medsmall)) ///
	ylab(`total_min'(`step')`total_max', labsize(medsmall) gmax) yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5) legend(off)
	gr export `graph_name'.eps, replace
}
end

******************************************************************
** Replication of Figure 5 (Arcenaux et al 2016) *********
** The Effect of Days to Election Conditional on Democratic Vote *
******************************************************************

tempfile twoway oneway
use "data/FoxNews_Master.dta", clear
gen dvprop=dv/100
gen daysdv=daystoelection*dvprop
gen days2dv=daystoelection2*dvprop
gen days3dv=daystoelection3*dvprop

global varlist "daystoelection daystoelection2 daystoelection3 dvprop daysdv days2dv days3dv Retirement seniorit qualchal_lag qualchal spendgap_lag spendgap distpart_lag RegPass Susp OtherPass Amend ProPart" 

logit PartyVote $varlist if PresencePartyUnity==1 & Republican==1 & FoxNews==1, cluster(dist2)

mat beta=e(b)[.,e(depvar) + ":"]
mat vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]

keep if e(sample)	
save current, replace										/* Save the sample */
 
drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace

cap mata: mata drop me()
mata
real matrix me(coef, x, x_new) {
	dydx=logistic(coef*x_new') - logistic(coef*x')			/* Replace logistic() with the appropriate link function as needed */
	return(dydx)
}
end

* Contour Plot
use current, clear
														
foreach x of varlist qualchal qualchal_lag Retirement {		/* Creating the necessary datasets */
	qui sum `x'
	replace `x'= (r(mean)>0.5)
}

foreach x of varlist seniorit spendgap_lag spendgap distpart_lag {
	qui sum `x'
	replace `x'=r(mean)
}

* Find the modal vote type
local dummies Amend OtherPass ProPart RegPass Susp
egen baseline = rowmax(`dummies')
replace baseline = 1-baseline
tabstat `dummies' baseline, save
mata props = st_matrix("r(StatTotal)")
mata st_local("modal", st_matrixcolstripe("r(StatTotal)")[selectindex(props :== max(props))[1,1],2])
foreach v in `dummies' {
	replace `v'=0
}
replace `modal'=1 if "`modal'"! = "baseline"

replace dvprop=round(dvprop,0.02) 
replace daystoelection=round(daystoelection,10)
collapse (mean) Amend OtherPass ProPart qualchal qualchal_lag RegPass Retirement Susp seniorit spendgap_lag spendgap distpart_lag (count) obs=Amend, by(daystoelection dvprop)
gen daystoelection2=daystoelection^2
gen daystoelection3=daystoelection^3
gen daysdv=daystoelection*dvprop
gen days2dv=daystoelection2*dvprop
gen days3dv=daystoelection3*dvprop

egen group_id=group(dvprop)
putmata wt=obs group_id=group_id X=($varlist 1), replace
preserve 
replace daystoelection = daystoelection+1
replace daystoelection2=daystoelection^2 
replace daystoelection3=daystoelection^3
replace daysdv=daystoelection*dvprop
replace days2dv=daystoelection2*dvprop
replace days3dv=daystoelection3*dvprop
putmata Z=($varlist 1), replace
restore
mata: mebr=me_byrow(coef, X, Z)
getmata (dydx_dx1 lb ub)=mebr

gen significant=(lb>0 & ub>0)|(lb<0 & ub<0)
qui sum daystoelection [fw=obs]
loc x1max: disp %9.1f r(max)
loc x1min: disp %9.1f r(min)
loc s1=round((`x1max'-`x1min')/4, 0.1)
qui sum dvprop [fw=obs]
loc x2max: disp %9.1f r(max)
loc x2min: disp %9.1f r(min)
loc s2=round((`x2max'-`x2min')/4, 0.1)

gen counter=_n
qui sum counter
loc extra1=`=r(max)'+1
loc extra2=`=r(max)'+2
loc extra3=`=r(max)'+3
loc extra4=`=r(max)'+4
set obs `extra4'

qui sum obs
replace significant=1 in `extra1'/`extra2'
replace significant=0 in `extra3'/`extra4'
replace obs=`=r(min)' in `extra1'/`extra4'
replace obs=`=r(max)' in `extra2'/`extra3'

qui sum dydx_dx1, detail
loc locut=`r(min)'+(`r(max)'-`r(min)')*1/5 
loc lmedcut=`r(min)'+(`r(max)'-`r(min)')*2/5
loc hmedcut=`r(min)'+(`r(max)'-`r(min)')*3/5
loc hicut=`r(min)'+(`r(max)'-`r(min)')*4/5

loc colr="navy*.5 ltblue*.5 white*.5 orange*.5 red*.5"

twoway (contour dydx_dx1 daystoelection dvprop if dydx_dx1!=., ccuts(`locut' `lmedcut' `hmedcut' `hicut') ccolors(`colr')) /*
*/ (scatter daystoelection dvprop [fw=obs] if significant==0, msymbol(oh) mlcolor(black%95) mlwidth(vthin) msize(*.25)) /*
*/ (scatter daystoelection dvprop [fw=obs] if significant==1, msymbol(o) mfcolor(black%95) mlwidth(none) msize(*.25)), /*
*/ xtitle("Democratic Vote Share", size(medsmall)) ytitle("Days to Election", size(medsmall)) ztitle("") ylab(`x1min'(`s1')`x1max', labsize(medsmall) grid gmax) xlab(`x2min'(`s2')`x2max', labsize(medsmall) grid gmax) xsize(7) ysize(5.5) /*
*/ zlabel(minmax `locut' `lmedcut' `hmedcut' `hicut', format(%9.6f) labsize(medsmall)) legend(off) clegend(title("Effect Size", size(medsmall) pos(12) justification(right)) width(5) height(25)) name(yx, replace) 
gr export Figure_5.eps, replace	 /* We do not use the contourp function here since the effects change signs and labels do not need to be rounded. */

******************************************************************
** Replication of Figures 7b and Appendices C1 and C2 ************
******************************************************************

** DAME with Bins by Deciles and Terciles, and Bins of Equal Width
use current, clear
gen obs=1
xtile group_id = dvprop, nq(10)
xtile group_id2 = dvprop, nq(3)
qui sum dvprop
loc mn=r(min)
loc mx=r(max)
loc atleast=100
tempvar size
gen group_id3=.
foreach j of numlist 20/3 {
	replace group_id3=ceil(`j'*(dvprop-`mn')/(`mx'-`mn'))
	egen `size'=total(obs),by(group_id3)
	qui sum `size'
	if r(min) >= `atleast' {
		continue, break
	}
	drop `size'
}

egen midpoint=median(dvprop),by(group_id)
egen midpoint2=median(dvprop),by(group_id2)
egen midpoint3=median(dvprop),by(group_id3)

putmata wt=obs group_id=midpoint group_id2=midpoint2 group_id3=midpoint3 X=($varlist 1), replace

preserve 
replace daystoelection=daystoelection+1
replace daystoelection2=daystoelection^2 
replace daystoelection3=daystoelection^3
replace daysdv=daystoelection*dvprop
replace days2dv=daystoelection2*dvprop
replace days3dv=daystoelection3*dvprop
putmata Z=($varlist 1), replace
restore
mata: dame=me_wt(coef, X, Z, group_id, wt)
mata: dame2=me_wt(coef, X, Z, group_id2, wt)
mata: dame3=me_wt(coef, X, Z, group_id3, wt)

* Mean case
use current, clear
qui sum dvprop
loc mn=r(min)
loc mx=r(max)

* Find the modal vote type
local dummies Amend OtherPass ProPart RegPass Susp
egen baseline = rowmax(`dummies')
replace baseline = 1-baseline
tabstat `dummies' baseline, save
mata props = st_matrix("r(StatTotal)")
mata st_local("modal", st_matrixcolstripe("r(StatTotal)")[selectindex(props :== max(props))[1,1],2])

collapse (mean) qualchal qualchal_lag Retirement daystoelection seniorit (median) spendgap_lag spendgap distpart_lag
foreach v in `dummies' {
	gen `v'=0
}
replace `modal'=1 if "`modal'"! = "baseline"

expand 21
gen dvprop=`mn' + (_n-1)*(`mx'-`mn')/20
gen daystoelection2=daystoelection^2 
gen daystoelection3=daystoelection^3
gen daysdv=daystoelection*dvprop
gen days2dv=daystoelection2*dvprop
gen days3dv=daystoelection3*dvprop
putmata X=($varlist 1), replace
preserve 
replace daystoelection = daystoelection+1
replace daystoelection2=daystoelection^2 
replace daystoelection3=daystoelection^3
replace daysdv=daystoelection*dvprop
replace days2dv=daystoelection2*dvprop
replace days3dv=daystoelection3*dvprop
putmata Z=($varlist 1), replace
restore
mata: mem=me_byrow(coef, X, Z)
getmata (dydx_mdx1 lbm ubm)=mem 
getmata (midpoint obs dydx_dx1 lb ub)=dame, force
getmata (midpoint2 obs2 dydx_dx2 lb2 ub2)=dame2, force
getmata (midpoint3 obs3 dydx_dx3 lb3 ub3)=dame3, force
save `oneway', replace
 
* Graphs
dame_plot `oneway' "DAME of Days to Election" "Democratic Vote Share" dvprop midpoint 5 Figure_7b

dame_plot `oneway' "DAME of Days to Election" "Democratic Vote Share" dvprop midpoint 5 Appendix_C2
gr_edit .title.text = {}
gr_edit .title.text.Arrpush Deciles
gr_edit .yaxis1.reset_rule -.0001 .0004 .0001 , tickset(major) ruletype(range)
gr rename FigureC12, replace

use `oneway', clear
replace dydx_dx1=dydx_dx2
replace lb=lb2
replace ub=ub2
replace obs=obs2
replace midpoint=midpoint2
save `oneway', replace
dame_plot `oneway' "DAME of Days to Election" "Democratic Vote Share" dvprop midpoint 5 Appendix_C1
gr_edit .title.text = {}
gr_edit .title.text.Arrpush Terciles
gr_edit .yaxis1.reset_rule -.0001 .0004 .0001 , tickset(major) ruletype(range)
gr rename FigureC11, replace

use `oneway', clear
replace dydx_dx1=dydx_dx3
replace lb=lb3
replace ub=ub3
replace obs=obs3
replace midpoint=midpoint3
save `oneway', replace
dame_plot `oneway' "DAME of Days to Election" "Democratic Vote Share" dvprop midpoint 5 Appendix_C3
gr_edit .title.text = {}
gr_edit .title.text.Arrpush "Bins of Equal Width"
gr_edit .yaxis1.reset_rule -.0001 .0004 .0001 , tickset(major) ruletype(range)
gr rename FigureC13, replace

grc1leg2 FigureC11 FigureC12 FigureC13, legendfrom(FigureC13) loff rows(1) ytol1title xtob1title ycommon
gr export "Appendix_C.eps", replace

******************************************************************
** Replication of Figures 1, 3, and 7a (Golder 2006) *************
** The Effect of Polarization Conditional on Threshold ***********
******************************************************************

tempfile twoway oneway
use "data/interaction3.dta", clear
global varlist polarization threshold polarization_threshold seatshare seatshare_2 incompatibility asymmetry asym_seat
xtprobit pec $varlist , re i(ident) 

matrix beta=e(b)[.,e(depvar) + ":"]
matrix vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]
preserve
keep if e(sample)
save current, replace /* save the sample */
drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace
restore

cap mata: mata drop me()
mata
real matrix me(coef, x, z) {
	int_coef_names = ("polarization","polarization_threshold")				/* The coefficients used to compute the derivative of the linear prediction */
	nam = st_matrixcolstripe("beta")
	k = J(cols(int_coef_names),1,.)
	for (j=1; j<=cols(int_coef_names); j++) {
		k[j] = selectindex(nam[.,2]:==int_coef_names[j])
	}
	dydx = (coef[.,k]*z'):*normalden(coef*x') 								/* Replace normalden() with the derivative of the inverse link function as needed */
	return(dydx)
}
end

keep if e(sample)															/* Creating the necessary datasets */
foreach var of varlist seatshare incompatibility asymmetry {
	qui sum `var'
	replace `var'=`=r(mean)'
}

collapse (mean) seatshare incompatibility asymmetry (count) obs=seatshare, by(polarization threshold)

gen polarization_threshold=polarization*threshold
gen seatshare_2=seatshare^2
gen asym_seat=seatshare*asymmetry

putmata wt=obs X=($varlist 1) Z=(1 threshold), replace
mata: mebr=me_byrow(coef, X, Z)
getmata (dydx_dx1 lb ub)=mebr
save `twoway', replace

use current, clear
gen obs=1
qui sum threshold
loc mn=r(min)
loc mx=r(max)
xtile group_id = threshold, nq(10)
egen midpoint=median(threshold),by(group_id)
putmata wt=obs group_id=midpoint X=($varlist 1) Z=(1 threshold), replace
mata: dame=me_wt(coef, X, Z, group_id, wt)

* Mean case
collapse (mean) $varlist [fw=obs]
expand 21
replace threshold=`mn' + (_n-1)*(`mx'-`mn')/20
replace polarization_threshold=polarization*threshold
replace seatshare_2=seatshare^2
replace asym_seat=seatshare*asymmetry

putmata X=($varlist 1) Z=(1 threshold), replace
mata: mem=me_byrow(coef, X, Z)
getmata (dydx_mdx1 lbm ubm)=mem
mata: group = dame[.,1]
getmata (midpoint obs dydx_dx1 lb ub)=dame, force
save `oneway', replace

* Graphs
heatm `twoway' "Polarization" "Effective Electoral Threshold" polarization threshold 1 4 Figure_3a
contourp `twoway' "Polarization" "Effective Electoral Threshold" polarization threshold 1 4 Figure_3b
dame_plot `oneway' "DAME of Polarization" "Effective Electoral Threshold" threshold midpoint 4 Figure_7a

******************************************************************
** Replication of Figure 1 ***************************************
******************************************************************

use `oneway', clear
append using current, keep(threshold)
twoway (hist threshold if dydx_mdx1==., percent yaxis(2) fcolor(gs8%75) lwidth(vthin) lcolor(white%25)) ///
(line dydx_mdx1 threshold, lpattern(solid) lwidth(vthin) lcolor(black%95)) ///
(rline lbm ubm threshold, lpattern(dash) lwidth(vthin) lcolor(black%95)), /// 
xlab(0(10)40, labsize(medsmall)) ylab(-0.001(.003).008, labsize(medsmall) axis(1)) ylab(#4, labsize(medsmall) axis(2)) ysca(alt) ysca(alt axis(2)) legend(off) ///
xtitle("Effective Electoral Threshold", size(medsmall) height(4)) ytitle("ME of Polarization at Its Mean", size(medsmall) height(4)) ytitle("Percent", axis(2) height(4) size(medsmall)) ///
yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5)
gr export Figure_1.eps, replace

******************************************************************
** Replication of Figures 4 and 8 (Nagler 1991) ******************
** The Effect of Registration Rules Conditional on Education *****
******************************************************************

tempfile twoway oneway
use "data/scobit.dta", clear
drop if newvote==-1
global varlist closing neweduc educ2 cloeduc cloeduc2 age age2 south gov
probit newvote $varlist

mat beta=e(b)[.,e(depvar) + ":"]
mat vcov=e(V)[e(depvar) + ":",e(depvar) + ":"] 
keep if e(sample)
save current, replace 													/* Save the sample */

drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace 

cap mata: mata drop me()
mata
real matrix me(coef, x, x_new) {
	dydx=normal(coef*x_new')-normal(coef*x')							/* Replace normal() with the appropriate function as needed */
	return(dydx)
}
end
 
use current, clear
qui sum age
replace age=r(mean)
foreach v in gov south {
	qui sum `v',detail
	replace `v'=r(p50)
}
 
collapse (mean) age south gov (count) obs=vote, by(closing neweduc)
gen age2=age^2
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*neweduc^2
putmata X=($varlist 1), replace
preserve
replace closing = closing+1
replace cloeduc = closing*neweduc
replace cloeduc2 = closing*educ2
putmata Z=($varlist 1), replace
restore
mata: mebr=me_byrow(coef, X, Z)
getmata (dydx_dx1 lb ub)=mebr
save `twoway', replace

* DAME
use current, clear
collapse (count) obs=vote, by(closing neweduc age south gov)
gen age2=age^2
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*neweduc^2
putmata wt=obs group_id=neweduc X=($varlist 1), replace 
replace closing = closing+1
replace cloeduc = closing*neweduc
replace cloeduc2 = closing*educ2
putmata Z=($varlist 1), replace
mata: dame=me_wt(coef, X, Z, group_id, wt)

* Mean case
use current, clear
qui sum neweduc
loc mn=r(min)
loc mx=r(max)
collapse (mean) age closing (median) gov south
expand 21
gen neweduc=`mn' + (_n-1)*(`mx'-`mn')/20
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*educ2
gen age2=age^2

putmata X=($varlist 1), replace
preserve
replace closing = closing+1
replace cloeduc = closing*neweduc
replace cloeduc2 = closing*educ2
putmata Z=($varlist 1), replace
restore
mata: mem=me_byrow(coef, X, Z)
getmata (dydx_mdx1 lbm ubm)=mem
getmata (midpoint obs dydx_dx1 lb ub)=dame, force
save `oneway', replace

* Graphs
contourp `twoway' "Closing Date" "Education" closing neweduc 0 3 Figure_4b
dame_plot `oneway' "DAME of Closing Date" "Education" neweduc midpoint 3 Figure_8a

use `oneway', clear
append using current,keep(neweduc)
twoway (hist neweduc if mi(dydx_mdx1), percent disc yaxis(2) fcolor(gs8%75) lwidth(vthin) lcolor(white%25)) ///
(line dydx_mdx1 neweduc, lpattern(solid) lwidth(vthin) lcolor(black%95)) ///
(rline lbm ubm neweduc, lpattern(dash) lwidth(vthin) lcolor(black%95)), /// 
xlab(1(1)8, labsize(medsmall)) ylab(-.004(.001).001, labsize(medsmall) axis(1)) ylab(#6, labsize(medsmall) axis(2)) ysca(alt) ysca(alt axis(2)) legend(off) ///
xtitle("Education", size(medsmall) height(4)) ytitle("ME of Closing at Its Mean", size(medsmall) height(4)) ytitle("Percent", axis(2) height(4) size(medsmall)) ///
yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5)
gr export Figure_4a.eps, replace

******************************************************************
** The Effect of Education Conditional on Registration Rules *****
******************************************************************

tempfile twoway oneway
use "data/scobit.dta", clear
drop if newvote==-1
global varlist closing neweduc educ2 cloeduc cloeduc2 age age2 south gov
probit newvote $varlist

mat beta=e(b)[.,e(depvar) + ":"]
mat vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]
keep if e(sample) 
save current, replace

drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace

* Creating the necessary datasets
use current, clear
qui sum age
replace age=r(mean)
foreach v in gov south {
	qui sum `v',detail
	replace `v'=r(p50)
}
 
collapse (mean) age south gov (count) obs=vote, by(closing neweduc)

gen age2=age^2
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*neweduc^2
putmata X=($varlist 1), replace
preserve
replace neweduc = neweduc+1
replace educ2 = neweduc^2
replace cloeduc = closing*neweduc
replace cloeduc2=closing*neweduc^2
putmata Z=($varlist 1), replace
restore
mata: mebr=me_byrow(coef, X, Z)
getmata (dydx_dx1 lb ub)=mebr
save `twoway', replace
 
* DAME
use current, clear
xtile group_id = closing, nq(10)
egen midpoint=median(closing),by(group_id)
collapse (count) obs=vote, by(closing neweduc age south gov midpoint)
gen age2=age^2
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*neweduc^2
putmata wt=obs group_id=midpoint X=($varlist 1), replace
replace neweduc = neweduc+1
replace educ2 = neweduc^2
replace cloeduc = closing*neweduc
replace cloeduc2=closing*neweduc^2
putmata Z=($varlist 1), replace
mata: dame=me_wt(coef, X, Z, group_id, wt)

* Mean case
use current, clear
qui sum closing
loc mn=r(min)
loc mx=r(max)
collapse (mean) age neweduc (median) gov south
expand 21
gen closing=`mn' + (_n-1)*(`mx'-`mn')/20
gen age2=age^2
gen educ2=neweduc^2
gen cloeduc=closing*neweduc
gen cloeduc2=closing*neweduc^2

putmata X=($varlist 1), replace
preserve
replace neweduc = neweduc+1
replace educ2 = neweduc^2
replace cloeduc = closing*neweduc
replace cloeduc2=closing*neweduc^2
putmata Z=($varlist 1), replace
restore
mata: mem=me_byrow(coef, X, Z)
getmata (dydx_mdx1 lbm ubm)=mem
getmata (midpoint obs dydx_dx1 lb ub)=dame, force
save `oneway', replace

* Graphs
contourp `twoway' "Education" "Closing Date" neweduc closing 0 3 Figure_4d
dame_plot `oneway' "DAME of Education" "Closing Date" closing midpoint 3 Figure_8b

use `oneway', clear
append using current,keep(closing)
twoway (hist closing if mi(dydx_mdx1), percent disc yaxis(2) fcolor(gs8%75) lwidth(vthin) lcolor(white%25)) ///
(line dydx_mdx1 closing, lpattern(solid) lwidth(vthin) lcolor(black%95)) ///
(rline lbm ubm closing, lpattern(dash) lwidth(vthin) lcolor(black%95)), /// 
xlab(0(10)50, labsize(medsmall)) ylab(.07(.01).12, axis(1) labsize(medsmall)) ylab(0(10)50, labsize(medsmall) axis(2)) ysca(alt) ysca(alt axis(2)) legend(off) ///
xtitle("Closing Date", size(medsmall) height(4)) ytitle("ME of Education at Its Mean", size(medsmall) height(4)) ytitle("Percent", axis(2) height(4) size(medsmall)) ///
yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) xsize(5.5) ysize(5.5)
gr export Figure_4c.eps, replace

******************************************************************
*Replication of Figures 6 and 7c (Robertson and Teitelbaum 2011) *
******************************************************************

tempfile twoway oneway
use "data/Robertson Teitelbaum 2011.dta", clear

global varlist "l_dispute open_penn l_gdp_pc_penn gdp_grth inflation_1 urban xratchg l_pop time"
global varlist1 "l_l_flows l_polity2 l_demflows l_dispute open_penn l_gdp_pc_penn gdp_grth inflation_1 urban xratchg l_pop time"

tsset country year 
gen l_l_flows=L.l_flows
gen l_polity2=L.polity2
gen l_dispute=L.dispute
gen l_demflows=l_l_flows*l_polity2
 
xtnbreg dispute $varlist1, re

mat beta=e(b)[.,e(depvar) + ":"]
mat vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]

keep if e(sample)
save current, replace

drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace

cap mata: mata drop me()
mata 
real matrix me(coef, X, Z) {
	int_coef_names = ("l_l_flows","l_demflows") 	/* The coefficients used to compute the derivative of the linear prediction */
	nam=st_matrixcolstripe("beta")
	k=J(cols(int_coef_names),1,.)
	for (j=1; j<=cols(int_coef_names); j++) {
		k[j]=selectindex(nam[.,2]:==int_coef_names[j])
}
	dydx=(coef[.,k]*Z'):*exp(coef*X') 				/* Replace exp() with the derivative of the inverse link function as needed */
	return(dydx)
}
end

use current, clear 									/* Replace the values of the covariates with their in-sample means */
foreach y of global varlist {
	qui sum `y'
	replace `y'=r(mean)
}

* Round the values of the constitutive term to reduce overplotting
egen l_flows1=cut(l_l_flows), at(-9.22, -7.22, -5.22, -3.22, -1.22, 1.22, 3.22, 5.22, 7.22, 9.22)
egen flows_mean=mean(l_l_flows), by(l_flows1)
replace l_l_flows=flows_mean
collapse (mean) $varlist (count) obs=dispute, by(l_l_flows l_polity2) 
gen l_demflows = l_l_flows*l_polity2

putmata wt=obs X=($varlist1 1) Z=(1 l_polity2) , replace
mata: mebr=me_byrow(coef, X, Z)
getmata (dydx_dx1 lb ub)=mebr
save `twoway', replace

* DAME
use current, clear
xtile group_id = l_polity2, nq(4)
egen midpoint=median(l_polity2),by(group_id)
gen obs=1
putmata wt=obs group_id=midpoint X=($varlist1 1) Z=(1 l_polity2) , replace
mata: dame=me_wt(coef, X, Z, group_id, wt)
  
* Mean case
qui sum l_polity2
loc mn=r(min)
loc mx=r(max)
collapse (mean) $varlist1
expand 21
replace l_polity2=`mn' + (_n-1)*(`mx'-`mn')/20
replace l_demflows = l_polity2*l_l_flows

putmata X=($varlist1 1) Z=(1 l_polity2), replace
mata: mem=me_byrow(coef, X, Z)
getmata (dydx_mdx1 lbm ubm)=mem 
getmata (midpoint obs dydx_dx1 lb ub)=dame, force
save `oneway', replace

* Graphs
heatm `twoway' "ln(FDI Flows)" "Polity 2" l_l_flows l_polity2 0 2 Figure_6
dame_plot `oneway' "DAME of ln(FDI Flows)" "Polity 2" l_polity2 midpoint 3 Figure_7c


********************************************************
******************* Appendix E1 *****************
********************************************************
use "data/Robertson Teitelbaum 2011.dta", clear 
gen any_dispute = dispute > 0
gen polity2_l = L.polity2
gen growing = (l_flows+L.l_flows+L2.l_flows)>(L3.l_flows+L4.l_flows+L5.l_flows)
eststo lo: logit any_dispute i.growing##c.polity2_l

tabstat growing,by(polity2)

matrix beta=e(b)[.,e(depvar) + ":"]
matrix vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]
drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear
putmata coef=(*), replace

mata
polity2_l=range(-10,10,1)
nr=length(polity2_l)
X1=(J(nr,1,(0,1)),polity2_l,J(nr,1,0),polity2_l,J(nr,1,1))
X0=(J(nr,1,(1,0)),polity2_l,polity2_l,J(nr,1,(0,1)))
pr1=logistic(coef*X1')
pr0=logistic(coef*X0')
ef = pr1-pr0
e1 = mean(pr1)'
e0 = mean(pr0)'
means=mean(ef)'
ra=mm_quantile(ef, 1, (0.025 \ 0.975))'
res=(polity2_l,e0,e1,means,ra)
end
 
clear
getmata (z pr0 pr1 est lb ub)=res

twoway (line est z, lpattern(solid) lwidth(vthin) lcolor(black%95)) ///
(rline lb ub z, lpattern(dash) lwidth(vthin) lcolor(black%95)), /// 
ytitle("Effect of Increasing FDI on the Probability of Strikes", height(4) size(medsmall)) xtitle("Polity Score", height(4) size(medsmall)) ///
ylab(,labsize(medsmall)) xlab(,labsize(medsmall)) ///
yline(0, lwidth(vthin) lcolor(black%95) lpattern(solid)) legend(off) xsize(5.5) ysize(5.5)
gr export "Appendix_E1.eps", replace

timer off 1
timer off 2
timer list